function x = LSF_LM(J, r, x0, lb, ub, lambda, tol, max_iter)
% Solve optimization problem using Levenberg-Marquardt algorithm
%
% Input arguments:
% J        - The Jacobian matrix
% r        - The residual vector
% x0       - The initial guess for the solution
% lambda   - The initial damping parameter
% tol      - The tolerance for the relative change in the solution
% max_iter - The maximum number of iterations
%
% Output arguments:
% x - The solution vector

x = x0';
iter = 0;
J_nan = isnan(J);
if sum(J_nan(:)) ~= 0
    return;
end
while iter < max_iter 
    

    iter = iter + 1;
    
    % Compute the gradient and Hessian
    g = J'*r;
    H = J'*J;
    
    % Add damping parameter to the diagonal of the Hessian
    H_lm = H + lambda*diag(diag(H));
    
    % Solve for the search direction
    p = - H_lm \ g;
    
    % Compute the new solution
    x_new = x + p;

    x_new = max(min(x_new, ub), lb);
    
    r_new = r + J*p;
    
    % Compute the relative change in the solution
    rel_change = norm(x_new - x) / norm(x);
    
    % Update the solution and the damping parameter
    if rel_change < tol
        break;
    elseif norm(r_new) < norm(r)
        lambda = lambda / 10;
        x = x_new;
        r = r_new;
    else
        lambda = lambda * 10;
    end
end
x = x';
end

